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WORKSHOP ON THE CE/SE METHOD 


Abstract 

The space-time conservation element and solution element (CE/SE) method, which 
was originated and is continuously being developed at NASA Glenn Research Center, is a 
high-resolution, genuinely multidimensional and unstructured-mesh compatible numerical 
method for solving conservation laws. Since its inception in 1991, the CE/SE method 
has been used to obtain highly accurate numerical solutions for ID, 2D and 3D flow 
problems involving shocks, contact discontinuities, acoustic waves, vortices, shock/acoustic 
waves/vortices interactions, shock/boundary layers interactions and chemical reactions. 
Without the aid of preconditioning or other special techniques, it has been applied to both 
steady and unsteady flows with speeds ranging from Mach number = 0.00288 to 10. In 
addition, the method has unique features that allow for (i) the use of very simple non- 
reflecting boundary conditions, and (ii) a unified wall boundary treatment for viscous and 
inviscid flows. 

The CE/SE method was developed with the conviction that, with a solid foundation 
in physics, a robust, coherent and accurate numerical framework can be built without in- 
volving overly complex mathematics. As a result, the method was constructed using a set 
of design principles that facilitate simplicity, robustness and accuracy. The most important 
among them are: (i) enforcing both local and global flux conservation in space and time, 
with flux evaluation at an interface being an integral part of the solution procedure and 
requiring no interpolation or extrapolation; (ii) unifying space and time and treating them 
as a single entity; and (iii) requiring that a numerical scheme be built from a nondissi- 
pative core scheme such that the numerical dissipation can be effectively controlled and, 
as a result, will not overwhelm the physical dissipation. Part I of the workshop will be 
devoted to a discussion of these principles along with a description of how the ID, 2D 
and 3D CE/SE schemes are constructed. In Part II, various applications of the CE/SE 
method, particularly those involving chemical reactions and acoustics, will be presented. 
The workshop will be concluded with a sketch of the future research directions. 
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A History of CFD Development 


1. From 60’s to 70’s - many new methods were proposed and tested. Among them, the projection 
method was most successful. To date, most of commercial codes are based on the projection 
method and its branches, e.g., the SIMPLE family. 

2. From early 80’s to early 90’s - CFD methods for aerodynamics and shock capturing were under 
vigorous development. In particular, significant resources were devoted to developing upwind 
schemes. The cost-to- value ratio in developing upwind schemes, however, has been disappointing. 
In general, upwind schemes are not robust and require special treatments to overcome numerical 
difficulties when applied to flow problems with practical geometries. Extension to flows with 
complex physics requires modifications in flux function formulation, which in general involve ad- 
hoc assumptions and approximations. 3D solutions have been scant. By and large, the upwind 
schemes have not been adopted by the commercial CFD community. Further development of 
upwind methods has not been a focal point of CFD research activities since early 90’s. 

3. From early 90’s to now - The CFD issues have been shifted to (a) unstructured mesh method, 
(b) high-performance parallel computation, (c) high-accuracy unsteady flow simulation, e.g., 
computational aeroacoustics (CAA) and large eddy simulation (LES) and (d) extension to inter- 
disciplinary applications, e.g., combustion, plasma dynamics, flow with multiple phases, crystal 
growth, etc. It is clear that a successful CFD method must be able to address these issues. 

4. The CE/SE method is designed to meet the above needs and thus can provide the CFD com- 
munity a robust and efficient numerical framework for high fidelity flow simulation in practical 
applications. 



Two Basic Beliefs 


The present development is guided by the following two basic beliefs that set 
the present method apart from the established methods. 

• In order to capture physics more efficiently and realistically, the modeling focus 
must be placed on the original integral form of the space-time physical con- 
servation laws, rather than the differential form. The latter form follows from 
the integral form under the additional assumption that the physical solution is 
smooth, an assumption that is difficult to realize numerically in a region of rapid 
change, such as a shock. 

• With proper modeling of the integral and differential forms themselves, the 
resulting numerical solution should automatically be consistent with the prop- 
erties derived from the integral and differential forms, e.g., the jump conditions 
across a shock and the properties of characteristics. As a result, a numerical 
method can be greatly simplified if the derived properties are not explicitly used 
in numerical modeling. 



Physics Considerations 


Flux-conservation in space-time be enforced locally (i.e., down to a single cell) 
and globally (i.e., over the entire solution domain). This is a critical requirement 
for accurate flow simulations, particularly if they involve long marching times 
and/or regions of rapid change. 

Space and time be unified and treated as a single entity. 

A numerical scheme be built from a nondissipative core scheme such that the 
numerical dissipation can be effectively controlled and, as a result, will not 
overwhelm the physical dissipation. 

A multidimensional scheme be genuinely multidimensional, i.e., it is constructed 
without using the dimensional- splitting approach, such that multidimensional 
effects and source terms (which are not aligned with special directions) can be 
modeled more realistically. 



Mathematics Considerations 


For easy implementation and efficient computing, the method should use the 
simplest logical structures and approximation techniques whenever possible. 

Flux at an interface separating two conservation cells be evaluated internally 
without using any ad hoc flux model or any characteristics-based technique. 

To be compatible with unstructured meshes, the spatial meshes be built from 
triangles in the 2D case, and tetrahedrons in the 3D case. 

The 2D and 3D schemes should share with their ID versions the same design 
principles so that they are straightforward extensions of the ID versions and 
share with the ID versions virtually identical fundamental characteristics. 

Avoid the use of special techniques that may limit application of the method 
and even cause undesirable side effects. 



Sod’s Shock Tube Problems with Non-Reflecting Boundary Conditions 


Extrapolation Boundary Conditions: (i) 

W" = (wm)“ri /2 and (“m*)" = {umxYjZljl, m = 1 , 2 , 3; n = 1, 2 , 3, . . . 
if (j, n) is a mesh point on the right spatial boundary; and (ii) 

(«m)J = and ( u mx )J = (u mx ) , m = 1 , 2, 3; n = 1 , 2, 3, . . . 

if (j, n) is a mesh point on the left spatial boundary. 

Steady- State Boundary Conditions: 

(^ra)j' = (^ra)j and {Urax) j — •> ^ 1, 2, 3, . . . 

where (j, n) is any mesh point on the right or left boundary. 



c 

c 


Appendix A. A Sample Program 


implicit real*8 (a-h,o-z) 
parameter (nxd=1000) 

dimension q(3,nxd), qn(3,nxd), qx(3,nxd), qt(3,nxd), 
* s(3,nxd), vxl(3), vxr(3), xx(nxd) 

c 

c nx must be an odd integer, 

nx = 101 
it = 100 
dt = 0.4d-2 
dx = 0.1d-l 
ga = 1.4d0 
rhol = 1 . dO 
ul = 0 . do 
pi = 1 . do 
rhor = 0.125d0 
ur = O.dO 
pr = 0 . ldO 
ia = 1 
c 

nxl = nx + 1 
nx2 = nxl/ 2 
hdt = dt/2 . dO 
tt = hdt*df loat (it ) 
qdt = dt/4.d0 
hdx = dx/2.d0 
qdx = dx/ 4 . dO 
dtx = dt/dx 
al = ga - l.dO 
a2 = 3. dO - ga 
a3 = a2/2 . dO 
a4 = 1 . 5d0*al 
u21 = rhol*ul 

u3l = pl/al + 0 . 5d0*rhol*ul**2 
u2r = rhor*ur 

u3r = pr/al + 0 . 5d0 *rhor*ur**2 
do 5 j = 1 , nx2 
q ( 1 , j ) = rhol 
q (2 , j ) = u2l 
q (3, j ) = u3l 
q(l,nx2+j) = rhor 
q(2,nx2 + j) = u2r 
q(3,nx2 + j) = u3r 
do 5 i = 1,3 
qx ( i , j ) = O.dO 
qx(i,nx2+j) = O.dO 
5 continue 

c 

open (unit = 8 , f ile= ' f or008 ' ) 
write (8,10) tt,it,ia,nx 
write (8,20) dt,dx,ga 
write (8,30) rhol, ul, pi 
write (8,40) rhor,ur,pr 
c 

do 400 i = 1 , i t 
m = nx + i - (i/2) *2 
do 100 j = 1 , m 
w2 = q (2 , j ) /q ( 1 , j ) 
w3 = q(3 , j ) /q (1, j ) 
f 2 1 = -a3*w2**2 
f 22 = a2*w2 

f31 = al*w2**3 - ga*w2*w3 
f32 = ga*w3 - a4*w2**2 
f33 = ga*w2 
qt(l,j) = - qx ( 2 , j ) 



qt ( 2 , j ) = - ( f 21*qx ( 1 , j ) + f22*qx(2,j) + al*qx(3,j)) 
qt ( 3 , j ) = - { £31*qx ( 1 , j ) + f32*qx(2,j) + f33*qx(3,j)) 
s(l,j) = qdx*qx(l,j) + dtx*(q(2,j) + qdt*qt(2,j)) 
s (2 , j ) = qdx*qx ( 2 , j ) + dtx* ( f 2 1* (q ( 1 , j ) + qdt*qt(l,j)) + 

* f 22* (q (2 , j ) + qdt*qt (2 , j ) ) + al*(q(3,j) + qdt*qt (3 , j ) ) ) 
s(3,j) = qdx*qx ( 3 , j ) + dtx* ( f 3 1* {q ( 1 , j ) + qdt*qt(l,j)) + 

* £32 * (q (2 , j ) + qdt *qt (2 , j ) ) + f 3 3 * ( q ( 3 , j ) + qdt*qt (3 , j ) ) ) 

100 continue 

if (i . ne . ( i/2 ) *2 ) goto 150 
do 120 k = 1,3 
qx(k,nxl) = qx(k,nx) 
qn (k, 1) = q (k, 1) 
qn(k,nxl) = q(k,nx) 

120 continue 
150 jl = 1 - i + (i/2) *2 
mm = m - 1 

do 200 j = 1 , mm 

do 200 k = 1,3 

qn (k, j + j 1) = 0 . 5d0* (q (k, j ) + q(k,j+l) + s(k,j) - s(k,j+l)) 
vxl (k) = ( qn ( k , j + j 1 ) - q(k,j) - hdt*qt (k, j ) ) /hdx 
vxr(k) = (q(k,j+l) + hdt*qt (k, j +1) - qn (k, j + j 1 ) ) /hdx 

qx(k,j+jl) = (vxl (k) * (dabs (vxr (k) )) **ia + vxr (k) * (dabs (vxl (k) ) ) 

* **ia) /( (dabs (vxl (k) )) **ia + (dabs (vxr (k) )) **ia + l.d-60) 
200 continue 

m = nxl - i + (i/2) *2 
do 300 j = 1 , m 

do 300 k = 1,3 

q (k, j ) = qn (k, j ) 

300 continue 
400 continue 


c 

m = nxl -it + (it/2) *2 
mm = m - 1 

xx (1) = -0 . 5d0*dx*df loat (mm) 
do 500 j = l,mm 
xx(j+l) = xx(j) + dx 
500 continue 

do 600 j = 1 , m 
x = q(2, j)/q(l, j) 

y = al*(q(3,j) - 0 . 5d0*x**2*q ( 1 , j ) ) 

z = x/dsqrt (ga*y/q(l, j ) ) 

write (8,50) xx ( j ) , q (1 , j ) , x, y, z 


600 

continue 




C 

close (unit=8) 



10 

format ( ' 

t = ' 

, gl4 . 7 , ' it = ' , 

i4 , ' ia 

20 

format ( ' 

dt = 

' ,gl4 .7, ' dx = ' 

,gl4 .7, ' 

30 

format ( ' 

rhol 

= ' ,gl4 .7, ' ul = 

' ,gl4.7 

40 

format ( ' 

rhor 

= ' , gl4 . 7, ' ur = 
f 8 . 4 , ' rho = ' , f 8 

' ,gl4.7 

50 

format ( ' 

x 

. 4 , ' u = 


* > 

M =' , 

f8 .4) 



stop 

end 





gamma = 
pl = 

, ' pr = 

' ,f8.4, ' 


nx = ' , i4 ) 
' ,gl4.7) 
,gl4.7) 
,gi4 . 7) 
p =' , f 8 . 4 , 
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Figure 1. 


The CE/SE solution of the extended Sod’s problem 
using the boundary conditions Eqs. (2.68) and (2.69) 
(At = 0.004, Ax = 0.01, CFL a 0.88, e = 0.5, a = 1) 
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Figure 1 (Continued). 
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Figure 2. — The CE/SE solution of the extended Sod’s problem 
using the boundary condition Eq. (2.70) 

(At = 0.004, Ax = 0.01, CFL « 0.88, e = 0.5, a = 1) 






The a Scheme 


Consider the PDE 


(i.i) 


du du 

m +a di~° 


where a is a constant. Let = x and X 2 = t be the coordinates of a 2D Euclidean space Eo ■ Because 
Eq. (1.1) <£> 

d(au) du 
v ; +-- = 0 


dx dt 


( 1 . 2 ) 

Gauss’ divergence theorem => 

(1.3) 
where 

(1.4) /i = (au, u), and ds = do n 


/ 

Js(V) 


h • ds = 0 



(a) Here (i) S(V) is the boundary of a space-time region V in E-z] and (ii) da and n are the area and 

the outward unit normal of a surface element on S(V). 

— # — * 

(b) h • ds is the space-time flux of h leaving the region V through the surface element ds. Thus 
Eq. (1.3) simply states that the total flux leaving the space-time region V through its boundary 
is equal to zero. 




5t(o ,n; 






n-1/2 ) 


n> n.) 











For any (x,t) 6 SE (j,n), let 


(1.5) u(x,t) & u* (x,t ; j,n) and h(x,t) & h* (x,t ; j,n) 
with 

(1.6) u*(x,t;j,n) d = uj + (■ u x )J(x - Xj) + (u t )*(t - t n ) 

(1.7) h*(x,t ;j,n) d = (au*(x,t]j,n),u*(x,t\j,n)) 

Here uj , (u x )™, and (itf)™ are constants in SE(j, n). 

Let u = u* (x,t ; j, n) satisfy the convection equation. Then 
(1-8) (u t )] = -a(u x )] 

Eqs. (1.6) and (1.8) =4> 

(1.9) u*(x,t]j,n) = uj + (tix)j [(a; - £j) - a(t - t n )] , (x,t) 6 SE(j,n) 

i.e., there are two independent marching variables u 1 -' and (u x )™ for each (j, n). 



Let AB G SE(j,n) be a part of S(V). Then 


(1.7) 

and 


h*(x, t jj, n) = f (uu*(x, t n), 


(1.9) u*(x,t ;j, n) = u" + (u x )" [(z - z 7 ) - a(t - t n ) 


(1.10) 


where 


h* -ds = 


AB 


h*{x c ,t c ; j,™) • ft 


X 


(i) (x c ,t c ) are the coordinates of the midpoint of AB. 


(ii) n is the outward unit normal of V at AB. 


h n)) 


, (x,t) <E SE(j,n) 


i 


(iii) £ is the length of AB. 



For any mesh point (j, n), let 


(l.n) 


/ -j- \ def AX . v n , def 0>At 

{ u x)j = ~ -r\ u * )j and 1/ = 


AX 


Hereafter the superscript symbol is used to denote a normalized parameter. 
Using Eq. (1.11) and the assumption j_± _L 

rh ' 


( 1 . 12 ) 

one has 

(1.13) 
and 

(1.14) 


i 


h* -ds = 0 


S(C£±0») 


CE_t h n) 


IT 




[(i - i/)u + (i - i/ a )«+]* = [(i - 1/)« - (i - «'*)«i']” +1 1 / / 2 2 ^ 






n 


i 


' ft 


— ^”-1 


[(1 + u)u - (1 - u 2 ) u+]” = [(1 + v)u + (1 - »' 2 )u£]“_ 1 1 / / J 2 


for all (j, n). To simplify notation, in the above and hereafter we adopt a convention 
that can be explained using the expression on the left side of Eq. (1.13) as an example, 

[(1 - u)u + (1 - «/>+]* = (1 - u)uj + (1 - «/*)(u+)! 


)" 

h 



Let 


(1.15) 

(«+)"+i/2 = f [“ - (! + 

(1.16) 


(1.17) 

« + )? = l [(»+);;# - 


By adding Eqs. (1.13) and (1.14) together, and using the above definitions, one has 

l 1 - 18 ) u ] = \ [(! “ * / )( s +)"+i 1 /2 +( 1 + v )( s ~) n j-ill 

Let 1 — v 2 / 0, i.e., 1 — v ^ 0 and 1 + i/ ^ 0. Then Eqs. (1.13) and (1.14) can be divided by 
( 1 — /') and (1 + n), respectively. By subtracting the resulting equations from each other, one has 

(1.13) («+)? = «)? 

The a scheme is formed by Eqs. (1.18) and (1.19). 

Note that the superscript symbol “a” that appears in Eqs. (1.17) and (1.19) is introduced to 
indicate that Eq. (1.19) is valid for the a scheme. 



(1.18) u? = \ [(1 - *)(•+)£$ + (1 + 

(1.19) («+)? = « + )j = i [(«+);;$ - (*-);:$' 

(1.9) u*(x,t\j,n) = Uj + (u x )J [(x - Xj) - a{t- t n )], (x,t) € SE(j,n) 


i 




Special Features of the a Scheme 


• It has the simplest stencil. 


• Its two amplification factors are identical to those of the Leapfrog scheme. 
Therefore it is nondissipative within its stability domain: o < M < i. 

• It is a two-way marching scheme, i.e., the forward marching scheme can be 
inverted to become the backward marching scheme and vice versa. 

• For each mesh point ( j,n ), it is associated with two discrete equations and two 
independent unknowns u J and {u x )J. As a result, the a scheme is consistent 
with a system of two PDEs involving two dependent variables with one of the 
PDEs being Eq. (1.1). 




In its construction, the surface integration over any interface separating two 
neighboring CEs is evaluated using the information from a single SE (no ex- 
trapolation or interpolation is used!). As a result, the local conservation rela- 
tion Eq. (1.12) leads to a global flux conservation relation, i.e., the total flux 
of h* leaving the boundary of any space- time region that is the union of any 
combination of CEs will also vanish. ^ n 


B 

(j-i/Zj-n-i/z) 



F 


C dfi/2 y n-1/z ) 


c 


V 




The <2-/2 Scheme 


The a-fi scheme is a solver of 


( 1 . 20 ) 


du du d 2 u 

dt + a di ~ = ° 


where a and fi > 0 are constants. Note that the integral form of the above PDE is 


(1.3) 


/ 

JS(V) 


h • ds = 0 


with 


( 1 . 21 ) 


h = ( au — fidu/dx , u) 


The construction procedure of the a-fi scheme is identical to that of the a scheme 
except that Eq. (1.7) is replaced by 


(1.22) h*(x,t ]j,n) d = ( au*(x,t \j,n) — ixdu*(x,t;j,n)/dx, u*(x,t ; J,n)) 



The a-e Scheme 


For the a-e scheme, the less stringent conservation condition 
(1.23) <f h*-ds = 0 

JS(CE(j,n )) 

is imposed at each (j, n). Again the local conservation condition Eq. (1.23) leads to 
a global conservation condition, i.e., the total flux ofh * leaving the boundary of any 
space-time region that is the union of any combination of CE(j,n), (j, n) £ Q, will 
also vanish. 


Because Eq. (1.23) 


(1.18) 


u? = - 

3 O 


1 r 


(1 - l/)(s + )J +1 1 / / 2 2 + (1 + 


1/2 


the a-e scheme shares with the a scheme the same first (principal) constituent equa- 
tion. 





c e. (j>; 


cE + fj>; 


C£-( j,n; 






For any (j, n), let 

(1-24) «&i/2 = u'll'l + (a</2)( U( )^; 2 2 

Using ( u t )" = and ^ = aAt/Ax , Eq. (1.24) => 

(1.25) 


W j±l/2 — [ U ^1^2 



n 


w- 


Because w'.^ 1 y 2 is a first-order Taylor’s approximation of n at (j ± 1/2, n). Thus 


(1.26) 


7 / ,n —n ,n . /.,'n —n ,n 

/'„.c+\n drf U j+l/2 U j-l/2 _ AX / u j+l/2 U j-l/2 

1 * " 4 ~ 4 ^ ax 


is a central-difference approximation of du/dx at (j, n), normalized by the factor ax/4. The a-e 
scheme is formed by Eq. (1.18) and 


(1.27) 

where e is a real number. 


(«+)? = «+)" + 2 *« + -< + )? 


(a) The superscript symbol “c” is introduced to indicate the central-difference nature of (u£ + )” . 

(b) The a-e scheme is stable if 0 < e < 1 and 0 < \v\ < 1; 

(c) The numerical dissipation associated with the a-e scheme generally increases with e. 

(d) The type of Numerical dissipation introduced by the second term in Eq. (1.24) is effective in 
damping out numerical instabilities that arise from the smooth region of a solution. But it is less 
effective in suppressing numerical wiggles that often occur near a discontinuity. 



The a-e-a-p Scheme 


Let 


(1.28) 


(», C +\ n ^— ^ (n, tU „. n \ — 

\ U x+)j — 2 ' j + 1/2 ~ U j) ~ 


AX 


u 


I n 


J+l/2 


— U 


n 


ax/2 



(1.29) 


Kt)? = ?( 


def 


1 


li 


n 



— 7 l ,n 

U J U j- 1/2 


ax/2 


Then 

(1-30) « + )? = 1 [(«St)" + («£)?] 


Let the function Wo be defined by (i) W o (0,0, a) = 0 and (ii) 



(1.31) 


W 0 (x_,x + ; a) 


\x+\ a X- + \X-\ a X + 
|x+| a + |x_| a ’ 


(|x + | + |x_| > 0) 


where x+, x_ and a > 0 are real variables. Note that W 0 (x_, x+; a), a nonlinear 
weighted average of x_ and x + , becomes their simple average if a = 0 or |x_| = |x+|. 



The a-€-<y-/3 scheme is formed by Eq. (1.18) and 


(1.32) 


Here 


(1.33) 


(«?)“ = K + Yj + M'C - K + y; + n« + - 


l«+ 


<•+ 


u'j+v 1 




«+)" =W C («+)?, «i)» 


Note that: 

(a) ( u "' + ) ; " = «+)" if a = 0 or (n c + + )J = (u c x ±)J. 

(b) The a-f-a-fi scheme generally is stable if0<e<l,/3>0 and a > 0. 

(c) We have 


(1.34) 


+ \n _ 


(<)■ 


f« + )y 

« + )j ife 

[(u w+ )J if 6 


/i = 0 
1/2 and (3 
1/2 and /3 


0 

1 


(a) In the smooth region of the solution, (u'|)" and (u£t.)” are more or less equal and, as a result, 
(u w+ — u£ + )™ becomes negligible. 

(b) Generally, with the choice of a = 1 or a = 2, the numerical dissipation introduced is sufficient 
to suppress numerical wiggles. 
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(a). — The dual space-time mesh 


Figure 22. — Concept of dual space-time meshes 



A rectangular space-time region 
shared by CE_(l/2,l/2) and CE + (0,l/2) 

Figure 22. (concluded) 







The Euler a Scheme 


Consider a dimensionless form of the 1-D unsteady Euler equations of a perfect gas. Let /?, v, p, 
and ~ be the mass density, velocity, static pressure, and constant specific heat ratio, respectively. Let 

(2.1) ui = P, u 2 = pv , u 3 =p/(7 - 1) 4- (l/2)/w 2 , 


(2.2) 

fl = «2, 

(2.3) 

h = (7 - l) w 3 + (l/2)(3 - t)(u 2 ) 2 /ui , 

(2.4) 

/ 3 = r {u 2 u 3 /ui - (l/2)(7 - 1 )(u 2 ) 3 /(wi) 2 

Then the Euler equations 

can be expressed as 

(2.5) 

5/ttx i n Q 

~ =0, m = 1,2,3. 

ot ox 

The integral form of Eq. (2.5) in space-time E 2 is 

(2.6) 

(p h Tn • ds = 0, m — 1, 2, 3, 

Js(V) 


where h w — ( / m , u m ), rn == 1,2,3, are the space-time mass, momentum, and energy current density 
vectors, respectively. 



Let (i) 


(2.7) /m,fc = dfm/diik, m,k = 1, 2, 3 

and (ii) (u m )j be the numerical version of at any (j,n). 

Because f m and / m? £ are functions of u m , for any (j, n), we can and will define (fm)™ 
and to be the values of f m and f m ,k, respectively, when u m , m = 1,2,3, 

respectively, assume the values of (u m )^, m = 1,2,3. Furthermore, because / m , 
m = 1,2,3, are homogeneous functions of degree 1 in the variables u m , m = 1,2,3, 
we have (see p. 11 of Methods of Mathematical Physics Vol. II, by R. Courant and 
D. Hilbert, interscience, 1962) 


( 2 . 8 ) 


3 


(M? = £(/m, *)>*)? 

fc = 1 


Nol e that Eq. (2.8) is not essential in the development of the ID CE/SE Euler solvers. 
However, in some instances, it is used to recast some equations into more convenient 
forms. 



For any (.r.f) £ SE(j, n), u m (x,£), f m (x,t) and h m (x,t ), respectively, are approximated by 


(2-9) <„(M ;j,n) = + (iw) "(* - *;) + (u m< );(t - t n ) 

(2.10) /;(M ;i, n) - (/ m ); + (U)J(X - *;) + (/ mf )J(t - t n ) 


( 2 . 11 ) 







Here we assume that 

3 

(2.12) (/„„)"=' and (/.»<)“ = 

fc=l 


3 


fc=l 


(2.13) (UmO” = -(/».)" 

Note that (i) the two equations in Eq. (2.12) are the numerical analogues of the analytical relations 

in dfm ^dfnduic dfm ^ df m du k 

( - 14) = and = 

fc=l fc=l 

respectively, and (ii) assuming Eqs. (2.9) and (2.10), Eq. (2.13) is equivalent to 

du* m (x,t\j,n) df^{x,t;j,n) _ 

at 


(2.13) 



Note that, by their definitions, 

(i) {fm)j and (/ m> *.)”, m = 1,2,3, are functions of (u m )?, m = 1,2,3, 

(ii) (fmx)ji m = 1,2,3, are functions of (u m )™ and (u mx )J : m = 1,2,3, and 

(iii) {fmt)], rn = 1,2,3, are functions of (u m )? and rn = 1,2,3. 

Furthermore, item (ii) and Eq. (2.13) imply that (tt m t)j, m — 1,2,3, are also 
functions of (w m )J and (u mx )J, m = 1,2,3. As a result, the coefficients on the right 
sides of Eqs. (2.9) and (2.10) are functions of (u m )^ and (u mx )™, m = 1,2,3. 

For all (j, n), we assume that 

(2.16) <b h* m -ds = 0, m = 1,2,3 

Js(CE±(j,n )) 

Using Eqs. (2.9)— (2.11), Eq. (2.16) implies that 


In l 71 - ( u \ n ~ 1 ! 2 4- — [(A !/ 2 , / \u 

\ a m)j \ u m)j±i/2 a \ a mx)j± i/ 2 ' \ u mx ) j 


(2.17) 


± f(/m)" 


AX L 


Tl — 1/2 / J. ^71 


j±l/2 


- (/m)l 


± 


(A /) 2 ' 


4ax 


(fmt)j±ii2 + (/mt) 


n 

3 


= 0. 



For each mesh point (j, n), let (i) 


(2-18) = m = l,2,3 

(ii) u n } and (u+)™, respectively, be the 3x1 column matrices formed by (it m )™ and 
(“«)”. ™ = 1 ) 2 , 3 , 

and 

(iii) (F + )’i be the 3x3 matrix formed by (At/Ax)(f m ,k) 1 j, m,k — 1,2,3. 

Then with the aid of Eqs. (2.8), (2.12) and (2.13), one can rewrite Eq. (2.17) as 
a pair of matrix equations, i.e., for any mesh point (j, n), 

(2.19) [(/ - F+)u + (7 - (F+) 2 ) 3+]* = [(7 - C+)u - (7 - (F+) 2 ) 
and 

(2.20) [(/ + F+)u - (7 - (F+) 2 ) 3+]; = [(7 + F+)u + (7 - (F+) 2 ) 3+];" 1 /’ 
where / is the 3x3 identity matrix. 



Note that the flux conservation conditions related to the a scheme, i.e. , 


(i.1.3) [(1 - v)u + (\ - v 2 )^]” = [(1 - u)u - (1 - i ^ 2 )uJ]J +1 1 / / 2 2 

and 

(1.14) {(1 + v)u - (1 - v 2 )u+] ” = [(1 + v)u + (1 - 

share the same algebraic structure with their Euler counterparts, i.e., 

(2.19) [ {I - + (I - (F + ) 2 ) ut} 7 ' = [(I-F + )u- (7-(F + ) 2 )uJ];;;/ 2 2 

and 

(2.20) [(/ + F+)«-(/-(F+) J )5+]"= [(/ + F+)u+(/-(F + ) 2 )«J]": i 1 / / 2 2 

In fact, the former pair will become the latter pair if the symbols 1, is, u and u+ are replaced by I , 
F + . u and uf, respectively. As a result, Eqs. (2.19) and (2.20) will be solved by a procedure similar 
to that used earlier to extract Eqs. (1.18) and (1.19) from Eqs. (1.13) and (1.14). 

However, because' 

(i) matrix multiplication is not commutative, 
and 

(ii) the matrix (E + )" is a function of (u m )”, m = 1,2,3, while v is a simple constant, 

as will be shown shortly, the algebraic structure of the solution to Eqs. (2.19) and (2.20) is more 
complex than that of Eqs. (1.13) and (1.14). 



The addition of Eqs. (2.19) and (2.20) 


(2.21) 

7 = \ { + N + F + )?~] 

where 


(2.22) 

(*+);+7/2 d = f + )^]"; 1 1 / 2 2 

and 


(2.23) 

( s "'-);: 1 1 / / 2 = [* + (i - f + )« 7];: 1 1 / / 2 2 


-1/2 / 


Note that Eqs. (2.21)— (2.23) are the Euler counterparts of 


(1.18) 

7 = 5 [(! - -)(^ + ); + 7/ 2 2 + (i + -)(.-« 

(1.15) 

(^);7i/ 2 2 d = i- - (i + -)-7i;7v 2 2 

and 


(1.16) 

(•-) nV a a = h + o -^];:^ 2 

respectively. 




To obtain the second part of the solution to 


(2.19) [(I - F+)S + (I - {F+)*) %]" = [(J-F+)u-(/-(F + ) 2 )u+]" +1 1 / / 2 2 

and 


(2.20) [(7+f+)«-(/-(f+) 2 )u+]; = [(/+F+)a+(/-(f+) 2 )u+];: 1 I / / 2 2 

existence of the inverses of both [J±(F + )]^ must be assumed. Let Eqs. (2.19) and (2.20) be multiplied 
from the left by 

[(/-F+)"]- 1 and [(/ + F+)"]- 1 

respectively. Let the resulting expressions be subtracted from each other. Then, because 
(2.24) [(/ - F + )(I + F + )}] = [(/ + F + )(I - F + )}] = [/ - (F+) 2 ]" 

one obtains 


(2.25) 


«)? = (ST)" 


where 

(2.26) («!+)" = \(3 + - S.)J 

(2.27) (S + )] d = f [(/ - F+)"]- 1 [(/ - F+)5 -(I- (F+) 2 ) ut^l 


(2.2S) 


(•?-)" = [(/ + F+)"]' 1 [(/ + F+)u+(/-(F + ) 2 )^]": i 1 / / 2 2 



Note that 


(2.25) (tZ+)J = K + )? 

(2.2G) (ul + )" d 4 f i(S + -5_)" 

( 2 . 27 ) ( S + )" d = f [(/- F +)"]- 1 [(/ - F+)u- (I-(F+f) «£]££ 

and 

( 2 . 2 5 ) ( S _)" ^ [(/ + F +)"]- 1 [(/ + F + )u + ( J -( F + )*) a +];- 1 / ^ 
respectively, are the Euler counterparts of 

(i.i 9) («+)? = («;+)? 


(1.17) 

k + )? = l [(.+)£$ - (-)”:# 

(1.15) 

(«+)£# = [« ■ - a + ■')«*]££ 

and 


(1.16) 




Let the marching variables at the (n — l/2)th time level be given. Then Uj can 
be evaluated using 

(2.21) «7 = \ { [(/ ■ - F + )s + ]T + l'll + [( J + F 

Because [I ± F + ]™ is a function of it follows that 

(2.27) (S + )J [(/ - F+)J] - 1 [(/ - F+) u -{I- (F+) 2 ) 

(2.28) (S_)? d = [(I+F + );]~ 1 [(I + F + )u+{I-(F + ) 2 )^: Ijl 

and 

(2.26) = |(S+ - 5-)" 

can also be evaluated. Thus Eq. (2.21) and 

(2.25) (fi+)? = K + 7 


form an explicit marching scheme — the Euler a scheme. 



An Existence Theorem 


Let 

(i) v and c be the fluid speed and sonic speed, respectively (u and c are known functions of u m , 
m — 1.2,3 [1]). 

(ii) Vj and c", respectively, denote the values of v and c when u m , m = 1,2,3, respectively, assume 
the values of (u m )j, m = 1,2,3. 

and (iii) 


(2.29) 


= -<="). (•*)? = 


n ( \Tl ^cf n 

Tx v >' 


Then it is shown in [2] that 


[(/-F+)"]- 1 and [(/ + F+)"]- 1 

exist if and only if 

(2.30) < = 1,2,3 

Obviously that Eq. (2.30) must be true if, for all (j,n), the local Courant number 

u j = f max{|(i/i)7|,|(i/ 2 )”U(i/3)j|} < 1 


(2.31) 



The Euler a-e Scheme 


For the Euler a-e scheme, the less stringent conservation conditions 


(2.32) 


i 


S{CE(j,n )) 


h* m ■ ds = 0, m — 1,2,3 


is imposed at each (j, n). It can be shown that Eq. (2.32) 4^ Eq. (2.21). 

For any (j,n), let 

(2.33) S-i 1/2 = ill'll + (a</ 2 )(“i)"±i/ 2 = (« - 2 ^)"± i/2 

where (u<)J±i /9 is the column matrix formed by (u m f)j±i /2 > m = 1,2,3. Also let 


(2.34) 


f Tl 

/ -*c+\n ^ J + 1/2 

V )j 


->7' n / 17' n _ ,7' n 

J- 1/2 Arc / U j + l/2 U j — 1/2 


Arc 


Then the Euler a-e scheme is formed by Eq. (2.21) and 

(2.35) (t£)? = (^ + )" + 2 e (5'+ - «“+)" 
where e is a real number. 

Note that it has been shown numerically that (i) the Euler a-e scheme generally is stable if 

(2.36) 0 < e < 1, and uj < 1 for all (j, n) 

and (ii) the numerical dissipation associated with the scheme increases as the value of e increases. 



The Euler a-e-a-/3 Scheme 


For any (j,n) and any ra = 1,2,3, let 
(2-37) «++)" = i(XJ" + 1/2 - («»)?) = ^ ( 

and 

(2.3S) «+-)" i((« m )7 - (Oy 1/2 ) = ^ ( 

with (w , „,)j± 1 / 2 an< i ( u m)j being the rath components of ttj^/ 2 and u”, respectively. Let be 

the /nth component of (u£ + )j. Then it can be shown that 

(2-39) « + . )? = 5 [(«£+)? + («£-)"] 

Let (u x + )J be the 3x1 column matrix formed by 

W 0 («t+)".(«*™ + . -)";«) , m = 1,2.3 

The Euler a-e-a-(5 scheme is formed by Eq. (2.21) and 

(<%)" = K + )j + 2e(«y - *?)] + /3(«x + - sy )" 


(Um)] 1/2 

ax/2 


(Oj + r/2 ~ KQj 
ax/2 


(2.40) 




Figure 5. 


A spatial domain formed from congruent triangles, 
showing the spatial projections of the mesh points 




Figure 6. — (a) The CEs associated with G' (b) the CEs associated with C" 
(c) The relative positions of the CEs of successive time steps 



t 



Figure 1 1 . — (a) Conservation elements CEO) (j, k, n + 1/2), € = 1 , 2, 3, 
and j, k, n = 0, ±1 , ±2, — . (b) Solution elements SEW Q> k, n + 1/2), 
j, k, n = 0, ±1 , ±2, •••). 


E-11063 Chang 9pt/1000R nm from E-9180flg5 





Figure 1 2. — (a) Conservation elements 0 . k. n + 1), ( = 1 , 2, 3, 

j, k, = 1/3, 1/3 ±1 , 1/3 ±2, — , and n = 0, ±1, ±2, (b) Solution 

elements SEP) Q, k, n + 1), j, k = 1/3, 1/3 ±1 , 1/3 ±2, -, and n = 0, 
± 1 , ± 2 , •••). 
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Figure 21. — Spatial projection of part of a 3D space-time mesh, 
showing the construction of a CE 
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